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■ Abstract 

(-H , The Multiquadric Radial Basis Function (MQ) Method is a meshless collocation method with global 

basis functions. It is known to have exponentional convergence for interpolation problems. We descretize 
nonlinear elliptic PDEs by the MQ method. This results in modest size systems of nonlinear algebraic 
equations which can be efficiently continued by standard continuation software such as AUTO and CON- 
TENT. Examples are given of detection of bifurcations in ID and 2D PDEs. These examples show high 
accuracy with small number of unknowns, as compared with known results from the literature. 

Keywords: Continuation, elliptic PDEs, bifurcation analysis, multiquadric radial basis function 
method. 
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^ ; 1 Introduction 

Nonlinear multidimensional elliptic partial differential equations (PDEs) are the basis for many scientific and 
engineering problems, such as pattern formation in biology, viscous fluid flow phenomena, chemical reactions, 
crystal growth processes, etc. In these problems it is crucial to understand the qualitative dependence of the 
solution on the problem parameters. 

During the past two decades the numerical continuation approach has become popular for qualitative 
study of solutions to nonlinear equations, see e.g. |33|, [jlO|, Q and references therein. Several software 
packages, such as AUT097 || and content pa], are currently available for bifurcation analysis of systems of 
nonlinear algebraic equations and ODEs, with only limited bifurcation analysis for ID elliptic PDEs. For 2D 
PDEs, we mention the software package pltmg |IJ that allows to solve a class of boundary value problems 
on regions in the plane, to continue the solution with respect to a parameter and even to compute limit and 
branching points. This software combines a sophisticated finite element discretization with advanced linear 
algebra techniques. Numerical continuation for ID and 2D elliptic PDEs is currently an active research 
area, see e.g. |gl), §7|, ||, §> §> S @, and @ Ch 10] for reaction diffusion equations; and f|, 
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for CFD. The typical approaches used are based on the finite element or finite difference discretization 
of the PDEs. They result in very large (thousands or tens of thousands for 2D problems) systems of 
nonlinear algebraic equations with sparse matrices. The continuation process is typically based on the 
predictor-corrector algorithms that require solving nonlinear systems by the Newton type method at each 
continuation step. For the bifurcation analysis during the continuation process, one usually needs to compute 
at least few eigenvalues of the Jacobian matrix at each continuation step. The methods currently used both 
for the continuation and the corresponding eigenvalue problems are variants of Krylov subspace methods 
and recursive projection (RPM). Solving the resulting linear system and the eigenvalue problem require 
sophisticated algorithms and considerable computer resources (CPU time, memory, disk space, etc.). 

In this paper we report results of numerical experiments with continuation and detection of bifurcations 
for ID and 2D elliptic PDEs discretized by the Multiquadric Radial Basis Function (MQ) method. The MQ 
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method was first introduced for solving PDEs in 1990 by Kansa p2[ , [ p3[ . It is a meshless collocation method 
with global basis functions which leads to finite dimensional problems with full matrices. It was shown to 
give very high accuracy with a relatively small number of unknowns (tens or hundreds for 2D problems) . The 
corresponding linear systems can be efficiently solved by direct methods. This opens a possibility for using 
standard continuation software, such as auto and CONTENT, designed for bifurcation analysis of modest 
size problems. We also note that the MQ method does not require predetermined location of the nodes as 
the spectral method does. 

In Section |^ we summarize previous results on solving PDEs by the MQ method and our experiments 
with an eigenvalue problem. 

In Section |3| we formulate an adaptation of the MQ method for the discretization of the parametrized 
elliptic PDEs. 

In Section |J we present results of our numerical experiments with continuation of solutions and detection 
of bifurcations for ID and 2D elliptic PDEs. 
In Section ai we discuss our results. 



2 Review of multiquadric method for elliptic PDEs 

2.1 Summary of previous results 

The concept of solving PDEs using the radial basis functions (RBF) was introduced by Kansa in 1990 [p2| , 
|p3f . He implemented this approach for the solution of hyperbolic, parabolic, and elliptic PDEs using the 
MQ RBFs proposed by Hardy (20) for interpolation of scattered data. 

There exists an infinite class of RBFs. A radial basis function, f(x), x £ K™, depends only upon the 
distances between the nodes. A MQ RBF is gj(x) — ((x — Xj) 2 + c 2 ) 1 / 2 , where Xj is a reference node and 
Cj is a shape parameter. In the comprehensive study by Franke |13| , it is shown that MQ RBFs have the 
excellent properties for the interpolation of 2D scattered data. Among studied RBFs still only the MQ RBFs 
are proven to have the exponential convergence for the function interpolation |2§[ ] , p9]| . 

The numerical experiments by Kansa p2]| , p3|, and Golberg and Chen |l5j show high efficiency and 
very accurate solution with the MQ scheme. Kansa [^3| showed that MQ method yields a high accuracy 
for parabolic and elliptic PDEs. Example for the transient convection-diffusion problem with steep initial 
front demonstrated highly accurate solution by the MQ method with a small number of nodes even for large 
cell Peclet number Pe. Test cases with 20 nodes for the MQ method ran for diffusion coefficient D in the 
range from 10 _1 to 10 -3 . The corresponding cell Pe number was from 0.5 to 50.0. Exact and MQ solution 
are indistinguishable graphically (apparent difference less than 10~ 4 ) for D = 10 _1 and 10~ 2 , while small 
deviation (5%) was observed at D = 10~ 3 ,Pe = 50. No instability or wiggles was seen. Finite difference 
simulation with K = 200 nodes and optimal combination of the central and upwind differences for the case 
D = 10 1 , Pe = 5 resulted in the error of 3%, which was still several orders less accurate than the MQ 
method solution. 

In the numerical experiments with modeling the von Neumann blast wave Kansa |23| compared the exact 
solution and its derivatives with the MQ solution (35 nodes) and with finite difference ones (50, 500 and 
5000 nodes). The error in value and gradients of pressure, density and energy was 10 -6 or less for the MQ 
method, and in the range from 10~ 4 to 10~ 2 for the best finite difference result with 5000 nodes. 

Golberg and Chen [[16) showed that the solution of the 3D Poisson equation could be obtained with only 
60 randomly distributed nodes to the same degree of accuracy as a FEM solution with 71,000 linear elements. 

Sharan, Kansa, and Gupta |3(| showed that the MQ method yields very accurate solutions for elliptic 
PDE problems, including the biharmonic equation, and that the MQ approach is simple to implement on 
domains with irregular boundaries. Dubai et al. |ll[] noted many benefits of using MQ RBFs to solve the 
initial value problem for a 3D nonlinear equation for the collision of two black holes. The resulting discrete 
system has 2000 unknowns and was solved directly. 

Buhmann || showed that RBFs and, in particular, MQ RBFs are useful for constructing prewavelets and 
wavelets. Wavelets are most frequently used in time series analysis, but there are results for using wavelets 
to solve PDEs p2[ , |3C|| , As Buhmann points out, one can generate true wavelets by an orthonormalization 
process. The wavelets are an elegant way to achieve the same results as multi-grid schemes. The MQ 
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RBFs are attractive for prewavelet construction due to exceptional rates of convergence and their infinite 
differentiability. 

Franke and Schaback paper Q] provides the first theoretical foundation for solving PDEs by collocation 
using the RBF methods. 

Kansa and Hon |24| studied several methods to solve linear equations that arise from the MQ collocation 
problems. They studied the 2D Poisson equation, and showed that ill-conditioning of the system of equations 
could be circumvented by using the sub-structuring methods. 

Kansa |2j| introduced the concept of variable shape parameters Cj in the MQ scheme that appeared 
to work well in some cases. In the work by Kansa and Hon |24| , a recipe based upon the local radius of 
curvature of solution surface was found to perform better than a constant shape parameter MQ scheme. A 
simple variable shape parameter formula is based the local radius of curvature. Kansa and Hon[p4| tested 
the MQ method for the 2D Poisson equation with a set of exact solutions like F = exp(ax + by), cos(ax + 
by), sm(ax + by), \og(ax + by + c), exp(— a(x — 1/2)2 — b(y— 1/2)2) and arctan(ax-|-6y). They showed obtained 
a high accuracy (up to 10 -5 ) and a small residual norm (10 -4 ) on a modest node size set (121 nodes) while 
locally adapting the shape parameter Cj. 

Franke |lj| compared (global) RBF interpolation schemes against many popular compactly supported 
schemes such as finite difference method, and found that the global RBF schemes were superior on six 
criteria. 

Madych |27| showed theoretically that the MQ interpolation scheme converges faster as the constant MQ 
shape parameter becomes progressively larger. 

The multi-zone method of Wong et al. ]3S| ] is yet another alternative method for improving computational 
efficiency. This method is readily parallelizable, and the conditioning of the resulting matrices are much 
better. 

Hon and Mao [£l| showed that an adaptive algorithm that adjusted the nodes to follow the peak of the 
shock wave can produce extremely accurate results in ID Burgers equation with only 10 nodes, even for 
extremely steep shocks with Re = 10 4 . 



2.2 A simple eigenvalue problem. 

Accurate approximation of eigenvalue problems is essential for bifurcation analysis of PDEs. We have not 
found references in literature on the MQ-solution of eigenvalue problems. We therefore present here results 
for an eigenvalue problem for ID Laplace operator. For details on the MQ discretization see Section [|. This 
is a scalar problem 
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that has the exact solution: 



(A m , U m (x)) = ((7rm) 2 , sin^rnx)) , m = 1, 2, ... 

where (A m , U m (x)) is the m — th eigenpair of ([l]). Introduce the mesh x n = nh, n = 0, 1, 
and consider the standard second order finite difference (FDM) discretization of ([!]): 

Un+l - 2u„ + U„„i 
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The corresponding approximate eigenpairs are given by 
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m = 1, ...,N - 1. 



(1) 



,N,h= l/N, 
(2) 



We also solved (|l|) using the MQ discretization for several values of the number K of internal nodes. Denote 

by 

( ^m^j ^rf^l ? ?7i — l, K the corresponding approximate eigenpairs. 
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Table 1: Eigenvalue problem: comparison of eigenvalues and eigenf unctions 
a) MQ method with uniform node distribution for K — 5, 7 and 9. 



TO 


A m (exact) 


^m^i K — 5 


Rel. crr.e* ly 


Rel.err e 1 ^ 


Rel. err.ej, if = 47 


1 


9.86961 


9.86596 


3.7 x 1(T 4 


3.7 x 10~ 4 


3.7 x 1CT 4 


2 


39.4784 


39.6492 


4.3 x 10-* 


5.2 x l(T ;l 


1.5 x 1CT 3 



777, 


A TO (exact) 


Kt W ,K = 7 


Rel. err. e^ uy 


Rel. err.e^ 


Rel. err.e^, K = 76 


1 


9.86961 


9.86821 


1.4 x 1(T 4 


9.9 x 10~ b 


1.4 x 10~ 4 


2 


39.4784 


39.4738 


1.2 x 10~ 4 


1.8 x 10~ 4 


5.7 x 10~ 4 


3 


88.8264 


89.3648 


6.0 x 10~ a 


1.1 x 1Q-' A 


1.3 x 10- 3 



m 


A m (exact) 




Rel. err.ef y 


Rel. err.e* ig 


Rcl. crr.e^, K = 117 


1 


9.86961 


9.86901 


6.0 x 10 _b 


5.0 x 10~ b 


6.0 x 10~ b 


2 


39.4784 


39.4846 


1.6 x 10~ 4 


2.1 x 10~ 4 


2.4 x 10~ 4 


3 


88.8264 


89.1667 


3.8 x 10~ a 


7.3 x 10" 3 


5.4 x 10~ 4 


4 


157.913 


159.689 


1.1 x 10~ 2 


2.5 x 10~ 2 


9.6 x 10~ 4 



b) MQ method with nonuniform node distribution for K = 7 and 9 



TO 


A m (exact) 


^ W ,K = 7 


Rel.err. e™ w 


Rel.err. % y 


Rel.err. e\, K = 3477 


1 


9.86961 


9.86961 


6.8 x 10~ 8 


3.0 x 10~ b 


6.8 x 10~ 8 


2 


39.4784 


39.4782 


3.2 x 10- B 


3.0 x 10~ 4 


2.7 x 10~ 7 


3 


88.8264 


88.8139 


1.4 x 10" 4 


6.5 x 10" 4 


5.4 x 10" 4 




TO 


A m (exact) 


\^,K = 9 


Rcl. err. e^ Jy 


Rel.err. e£/ y 


Rel.err. e\, K = 950 


1 


9.86961 


9.86960 


9.1 x 10~ Y 


2.3 x 10~ b 


9.1 x 10~ Y 


2 


39.4784 


39.4783 


1.4 x 10~ b 


2.0 x 10~ 5 


3.6 x 10~ b 


3 


88.8264 


88.8241 


2.6 x 10~ 5 


1.8 x 10~ 4 


8.2 x 10~ b 


4 


157.913 


157.882 


1.9 x 10~ 4 


1.8 x 10~ 3 


1.5 x 10~ 5 



The results of our computations are summarized in Table [l]. We use the notation e 1 ^ ® , for the relative 
errors in A^f 1 ^, A^, respectively, and the notation e£f for the Loo-norm error in C/jJf^ . For each MQ 
solution we provide a comparison with the FDM solution that has a sufficient number of nodes to give the 
same accuracy for Ai as the MQ method. In Part (a) of the table we use the uniform node distribution 
for the MQ method. Part (b) of the table shows that the accuracy of the MQ method can be significantly 
improved by adapting the node distribution: we moved only two nodes adjacent to boundary to reduce their 
distance from the boundary to hi = h/A (while the remaining nodes are distributed uniformly). 

One can see that the MQ method can give a highly accurate solution with a small number of unknowns, 
10 — 100 times smaller than the number of unknowns in the FDM for the same accuracy. 



3 Discretization of nonlinear elliptic PDEs by the MQ method 

We consider the second order system of n parametrized nonlinear elliptic partial differential equations 

D(a)Au- f(Vu,u,x,y,a) = 0, a £ R, u(-), /(•) 6 M™, (x, y) £ Q C M 2 , (3) 

where D{a) is a positive diagonal n x n matrix, that dependents smoothly on a, subject to boundary 
conditions 



0,/ h (-)e 



(4) 



</<> 



Here a is a control parameter, and we are interested in studying the dependence of the solutions to the 
boundary value problem (3), (4) on a. 
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We discretize the continuous problem by the multiquadric radial basis function (MQ) method , |23| , 
p8[ as follows. Introduce a set Qh of nodes (./V internal and Nt, on the boundary) 

0/j = {{xuVi) |i=i,jvC fi, (xi,y t ) \i=N+i,N+N b C- 
and look for the approximate solution to (ph, (4) in the form 

j=N-l j=N+N b 

u h (x,y)=a + ^2 a j (9j{cj,x,y) - g N (c N ,x,y)) + aj {gj(cj,x,y) - g N (c N ,x,y)) , (5) 

3=1 3=N+1 

where aj 6 K™ are the unknown expansion coefficients and 

gj(cj,x, y) = \J{x- x 3 ) 2 + (y- %) 2 + c 2 , j = 1, AT + 7V 6 , 

are the MQ basis functions, and Cj > is called a shape parameter [23]. We then substitute Uh{x,y) into 
(|J), (4) and use collocation at the nodes 8/,, to obtain a finite dimensional system 

<Pi(a,a) = D(a)Auh(xi,yi) - f(Vu h (xi,yi),Uh(xi,yi),Xi,yi,a) = G, i = 1, ...,iV, (6) 

^(a.a) = / 6 ( a " h ^' yi) , Mh (a; i , 2 / < ),a:i,yi,a) = 0, i = JV + 1, JV + N b . (7) 

We next modify the discretized system to make it more suitable for continuation and bifurcation analysis. 
1) We eliminate aj, j = N + 1, N + N^, associated with the boundary nodes, so as to minimize the number 
of unknowns. 2) We reformulate (5) in terms of nodal values Uj so that to have the correct eigenvalue problem 
(to avoid dealing with matrix stencils) for the Jacobian matrix of (Q) for detecting bifurcations during the 
continuation process. 

This is accomplished as follows. Denote a 1 — (ao, a\, ajv-i) £ K" xAr , a 2 = (ajy+i, ajv+jv b ) £ M. nxNb , 
(f = (<p 1: ...,<p N ), (fi b = (ip\, (fi b Nb ), and rewrite the system (|), (7) as 

(p{a\a 2 ) a)=0, ^)eM" xAr , (8) 

<p b (a\a 2 ,a) = 0, /(•) e R nxN K (9) 

Assuming that the implicit function theorem is applicable here (which is usually the case), we solve (9) for 
a 2 to obtain 

a 2 = %j){a , a), or, in components, aj = ip -,(a , a), j = N + 1, N + Nb- (10) 
Substituting this into (8) yields 

tp{a 1 ,^{a 1 ,a),a)=Q, (p{-)eR nxN . (11) 

We next want to reformulate (11) in terms of the nodal values U = (ui, U2, mat) G R™ x N of the approximate 
solution at the internal nodes defined by Ui = Uh(xi, t/,). To this end we first eliminate a 2 = (ajv+i, dN+N b ) 
from (||) by substituting (|l0|) into (Q) to obtain 

j=N-l j=N+N b 

u h (x,y)=a a + a i (9j{cj,x,y) - g N (c N ,x,y)) + ^ ip j (a 1 ,a) (gj(cj,x,y) - g N (c N ,x,y)) (12) 

3 = 1 3=N+l 

We now define the map V : a 1 \ — > U = T(a 1 ). For i = 0, ...,N - 1 : 

j=N-l j=N+N b 

u t =a + ^2 (gj( c j> x t,yt) ~ 9N(cN,x l ,y l ))a j + ^ (9j(cj, x t , y t ) - g N (c N , x t , y t )) ^-(a 1 , a),. (13) 

3=1 3=N+1 

Finally, substituting a 1 = T' 1 ^) into (11), we arrive at the finite dimensional continuation problem 

G(U,a)=0, U,G(-) G R nxN ,a 6 M, (14) 

where 

G(C/, a) = if (r- x (E/), V (r- x (C/), a) , a) , r : l w -> 1 N , € R nXiVi> - 
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Table 2: ID Gelfand-Bratu equation: limit point comparison 
a) results for MQ correspond to uniform node distribution 





|J, exact 


[ j, K = 800 


MQ(u), K = 5 


MQ(u), K = 7 


MQ(u), K = 9 


A 


3.513831 


3.5137 


3.512609 


3.514224 


3.514047 


rel. error 




3.7 x 10~ b 


3.5 x 10~ 4 


-1.1 x 10~ 4 


-6.1 x 10~ b 



b) results for MQ correspond to nonuniform node distribution 





[ 25 1 , K = 50 


[25|, K = 500 


MQ(nu), K = 5 


MQ(nu), K = 7 


MQ(nu), K = 9 


A 


3.51145 


3.51380 


3.514010 


3.513809 


3.513828 


rel. error 


6.8 x 10~ 4 


8.8 x 10~ b 


-5.1 x 10~ b 


6.3 x 10~ b 


8.5 x 10~ Y 



Remark 1 Note that in the case that the boundary condition (4) is linear, ipj are linear, and consequently 
r is an N x N matrix. 

In Section [| we consider examples of continuation of ID PDEs with fl = (0, 1) and 2D PDEs with 
fl = (0, 1) x (0, 1). In all 2D examples we have the same number N s of nodes in x and y directions. We 
choose a constant shape parameter Cj = s/(N s — 1). Our typical choice for s is 4 < s < 12. 

We use two types of node distributions. In the case of uniform node distribution (xk, yi) = (kh, Ih), 
k, I = 0, ...,N S , h = In the case of nonuniform node distribution, the nodes adjacent to the boundary 

dVl are placed at the distance h = hih from dfl, 0.1 < hi < 0.5, while the remaining nodes are distributed 
uniformly. A criteria for the choice of hi was a minimum of L2-norm of the residual in Q. 



4 Numerical experiments for 1-D and 2-D elliptic PDEs 

We present several examples of continuation of solutions to systems of nonlinear ID and 2D elliptic PDEs. 
Each problem is discretized by the MQ method described in Section ||. We then perform continuation of the 
resulting system of algebraic equations ([14]) with AUT097. The principal goal of our examples is to assess the 
accuracy of the detection of bifurcation points. We compare our results with some published results and, in 
one case, the results of our computations with an example in AUT097 and content. We will use throughout 
the notation K for the number of unknowns in a particular method. For our MQ method K = n x N , where 
n is the dimension of the system and N is the number of internal nodes. We denote by MQ(u) and MQ(nu) 
our MQ method with the uniform and nonuniform node distribution, respectively. 

Example 1 ID Gelfand-Bratu equation. This is a scalar problem 

u" + Ae" = 0, in ft = (0, 1), (15) 
u(0) = u(l) = 0, 

that appears in combustion theory and is used as the demo example exp in AUT097 (forth order adaptive 
orthogonal spline collocation method) and demo example in brg in CONTENT Jjjj/ (third order adaptive finite 
difference method). There is a limit (fold) point on the solution curve. We take the value of A at the limit 
point found from demo exp (K > 50) as exact. The following table ^ shows comparison between numerical 
results in H], our numerical results and our experiments with content. 



b (16) 



Example 2 ID Brusselator problem. This is a reaction diffusion model for a trimolecular chemical reaction 

jtu" - (b+ l)u + u 2 v + a = 0, + bu - u 2 v = 0, in fl =(0,1), 

u(0) = u(l) = a, v(0) = v(l) = i 

A stationary bifurcation occurs [6, Eq. (24)] at 

di 9 n 2 n 2 , I 2 a 2 
b n = 1 + -ra 2 + —p^d! + -^-r- > 0. 

02 l z -K z n z di 
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Table 3: ID Brusselator equation: bifurcation points comparison 



a) bifurcation point b\ 





exact 


HI, # = 80 


MQ(u), K=1Q 


MQ(u), K = U 


MQ(u), K = 18 


bi 


19.680174 


19.67547 


19.67366 


19.67786 


19.67919 


rel. error 




2.4 x 10" 4 


3.3 x 10" 4 


1.2 x 10" 4 


5.0 x 10" b 


b) bifurcation point b 2 




exact 


[6], K = 80 


MQ(u), K = 10 


MQ(u), K = 14 


MQ(u), A" = 18 


b 2 


48.681060 


48. 6004 


48.57476 


48.63168 


48.65605 


rel. error 




1.7 x 10~ a 


2.2 x IQ-' A 


1.0 x 10~ ;5 


5.1 x 10~ 4 



Table 4: ID pattern formation problem, bifurcation points 



|8j, numerical] 


0.047 


0.080 


0.093 


0.159 


0.140 


0.238 


0.186 


0.317 


0.232 


|8j, analytic] 


0.0465 


0.0793 


0.093 


0.159 


0.140 


0.238 


0.186 


0.317 


0.233 


MQ(nu) 


0.0465 


0.0793 


0.093 


0.159 


0.140 


0.238 


0.186 


0.317 


0.233 



For I = d\ = 1, d,2 = 2, a = 4, n = 1, 2 this gives simple bifurcations: bi = 9 + it 2 + \ = 19.680174, 
62 = 9 + 47r 2 + -^- = 48.681060, correspondingly. For the second order central difference method with uniform 
mesh of 41 mesh points (K — 80 unknowns), the corresponding approximate bifurcation points were found 
in Section 6.1]. The following table^ shows comparison between analytical, numerical results /|^, Section 
6.1] and our numerical results for values of b\ and b 2 at simple bifurcation points. 



Example 3 Pattern formation in a ID system with mixed boundary conditions 

fyu" + (3 - ku ~ uv 2 — 0, Sfyv" + nu + uv 2 - v = 0, inf!=(0,l) 

e 1 %± = p(i-o 1 )(e 3 u s -u), se 2 ^ = d P (i-e 2 )(e 3 v s -v), on on = {0,1}. [ ' 



Here 8i G [0,1], i — 1,2,3, are homotopy parameters. For d\ = 10~ 5 . ui = 10~ 2 . 8 = 0.14, (3 = 1.0, 
K = 0.001, (9i,02, O3) — (1,1,0) (Neumann problem). Eq. (17) was discretized by the second order central 
difference method with equidistant mesh of AX mesh points (K = 80 unknowns) . The following table Table 
1] shows a comparison between analytic and numerical results for values of I at simple bifurcation points. 



Our numerical results (MQ(nu) method) with K = 18, coincide with the analytic results above. In 
addition, we found a bifurcation point at I = 0.279. 



Example 4 2D Bratu problem 



Au + Xe ll ,n = (0,l)x(0,l), (18) 
u on = 0. 



This problem was studied in (35 J. It was discretized with the second order central difference method with 



uniform mesh and then continued using Implicit Block Elimination based on Recursive Projections. A limit 
point was detected for some value of X (not reported in the paper), and spurious limit points were detected 
for K = 961, 1521, 2209, 3025 and A sufficiently small. We reproduced the bifurcation diagram in no 
spurious limit points were detected. The following table || gives the values of A at the limit point computed 
by MQ method. 



Example 5 2D Brusselator problem. 

^■Au-{b+l)u + u 2 v + a = Q, ^Av + bu-u 2 v = 0, in O = (0, 1) x (0, 1), . , 

I I b \ ^) 

u \en= a, v \ 9a = -. 
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Table 5: 2D Bratu equat 


ion, limit point 




|35J, 225 < K < 3025 


MQ(u), K = 25 


MQ(u), K = 49 


MQ(u), K = 81 


A 


not reported 


6.873498 


6.840836 


6.827400 



Table 6: 2D Brussclator equation: bifurcation points, uniform node distribution for MQ 



a) bifurcation point b\ 





exact 


J5J, K= 800 


MQ(u), K = 50 


MQ(u), K = 72 


MQ(u), K = 98 




29.144494 


29.104774 


29.16280 


29.17050 


29.16062 


rel. error 




1.4 x 10" 3 


-6.3 x 10" 4 


-8.9 x 10~ 4 


-5.5 x 10" 4 


b) bifurcation point 6 2 




exact 


0, K= 800 


MQ(u), K = 50 


MQ(u), K = 72 


MQ(u), K = 98 




88.058156 


87.47325 


87.61578 


87.86924 


88.00143 


rel. error 




6.6 x IQ-' A 


5.0 x lO" 3 


2.1 x IQ~' 6 


6.4 x 10~ 4 



/ TO 2 


fnj 


a 2 


f P ) 


VP" 




7T 2 d2 


\m 2 + l 2 n 2 J 



A stationary bifurcation occurs Eq. (2.26)] fo 

6 m ,„ = 1 + ^-a 2 + d X Tx A ( '—r- + n z I + -j— [ ' I > 0. 

For I = di = 1, d,2 = 2, a = 4, (m,n) = (1,1), (m, n) = (2,2) £/us gives simple bifurcations: b\ t i — 
9 + 27r 2 + -ij, 62,2 = 9 + 87T 2 + 4j-, correspondingly. For the second order central difference method with 
equidistant mesh of 21 mesh points, the corresponding approximate bifurcation points are found in Section 
5]. The following tables show comparisons between analytical, numerical results Eq. (2.26)] and our 
numerical results for values 0/61,1 and 62,2 ok simple bifurcation points. 



A Hopf bifurcation occurs §, Eq. (2.26)] for 

b m ,n = 1 + a 2 + (di + d 2 ) + n 2 ^ n 2 

for some m, n, and I large enough. For / = 10, d% = d<z = 1, a = 10, (m,n) — (1,2), this gives a Hopf 
bifurcation at b 1>2 = 101 + 2 + 2 2 ) vr 2 = 180.15, see table |. 



5 Conclusions. 

We presented the results of our experiments with the MQ method for continuation of solution of nonlinear 
ID and 2D elliptic PDEs. We used small number of unknowns and obtained a high accuracy for detecting 
bifurcation points in our examples. Here are some sample results. 



Table 7: 2D Brusselator equation: bifurcation points, nonuniform node distribution for MQ 
a) bifurcation point b\ 





exact 


MQ(nu), K = 50 


MQ(nu), K = 72 


MQ(nu), K = 98 




29.144494 


29.14621 


29.14726 


29.14431 


rel. error 




-5.9 x 10" 5 


-9.5 x 10" 5 


6.3 x 10" B 


b) bifurcation point 62 




exact 


MQ(nu), K = 50 


MQ(nu), K = 72 


MQ(nu), K = 98 


62,2 


88.058156 


88.15470 


87.93391 


88.07288 


rel. error 




-1.1 x lO" 3 


1.4 x lO" 3 


-1.7 x 10~ 4 



8 



Tabic 8: 2D Brusselator equation, Hopf bifurcation point 





exact 


MQ(u), K = 50 


MQ(nu), K = 50 


MQ(u), K = 72 


MQ(u), K = 98 


61,2 


180. 15 


181.8625 


180.7880 


181.0696 


180.492 


rel. error 




-9.5 x lO" 3 


-3.5 x 10~ a 


-5.1 x 


-1.9 x lO" 3 



(i) For the limit point in the ID Gelfand-Bratu equation, the MQ method with 9 unknowns gives the 
relative errors 6.1 x 10 -5 and 8.5 x 10 for the uniform and nonuniform node distributions, respectively. 
The relative error in the finite difference method with 500 nodes is 8.8 x 10~ 6 . 

(ii) For the two bifurcation points in the 2D Brusselator problem, the MQ method with 98 unknowns 
gives the relative errors 5.5 x 10~ 4 , 6.4 x 10~ 4 for the uniform node distribution and 6.3 x 10 -6 , 1.7 x 10~ 5 
for the nonuniform node distribution. The corresponding relative errors in the finite difference method with 
800 nodes are 1.4 x 10" 3 , 6.6 x 10~ 3 . 

(iii) for the first in the eigenvalue problem for the ID Laplace operator with 9 unknowns gives the relative 
error 6 x 10~ 5 and 9 x 10~ 7 for the uniform and nonuniform node distributions, respectively. This is equivalent 
in accuracy to 117 and 950 node solution, respectively by the finite difference method. 

The increase of the number of unknowns results in a better accuracy but also in a larger condition number 
of the operator T mapping solution nodal values to the expansion coefficients. This condition number is a 
limiting factor in our experiments. In the future, we plan to implement the ideas of Kansa et al. |24| to 
circumvent this ill conditioning problem. 

In addition we found that even a simple adaptation of the nodes adjacent to the boundary can lead 
to a dramatic improvement of the accuracy in detecting bifurcation points. Adaptive choice of the shape 
parameter is another way to improve the accuracy that we plan to investigate. 

Our results show that MQ method is an efficient method for continuation of solution nonlinear PDEs. 

Acknowledgments. Support from the NASA grant NAG8-1229 is acknowledged by the first author. 
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